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PM ' Abstract 



The composition of ultra-high energy cosmic rays is an important issue in as- 
troparticle physics research, and additional experimental results are required 



O . 

Vh , for further progress. Here we investigate what can be learned from the sta- 

^ ' tistical correlation factor r between the depth of shower maximum and the 

muon shower size, when these observables are measured simultaneously for a 
set of air showers. The correlation factor r contains the lowest-order moment 

>• . of a two-dimensional distribution taking both observables into account, and 

it is independent of systematic uncertainties of the absolute scales of the two 



t^^ I observables. We find that, assuming realistic measurement uncertainties, the 

value of r can provide a measure of the spread of masses in the primary 

^ ■ beam. Particularly, one can differentiate between a well-mixed composition 

^ . (i.e., a beam that contains large fractions of both light and heavy primaries) 

and a relatively pure composition (i.e., a beam that contains species all of a 
similar mass). The number of events required for a statistically significant 

K> ' differentiation is ~ 200. This differentiation, though diluted, is maintained 

H ■ to a significant extent in the presence of uncertainties in the phenomenology 

■ ■ ■ of high energy hadronic interactions. Testing whether the beam is pure or 

well-mixed is well motivated by recent measurements of the depth of shower 
maximum. 
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1. Introduction 

The composition of ultra-high energy (UHE) cosmic rays is a key question 
in astroparticle physics research. There is relatively little guidance from 
theory about what particles to expect. Typically, particle masses ranging 
from protons to heavy nuclei, such as iron, are regarded. Experimental input 
is needed to clarify the situation. 

Measurements of the depth of shower maximum Xmax from the Pierre 
Auger Observatory [l] show a statistically significant flattening of the elon- 
gation rate near 2 x 10^* eV. In addition, the natural shower-to-shower fluc- 
tuations of Xmax are reported to decrease from approximately 60 g/cm^ at 
2 X 10^^ eV to 30 g/cm^ at 4 x 10^^ eV. These observations suggest an increase 
of the average nuclear mass of UHE cosmic rays with energy. In particular, 
the shower-to-shower fluctuations observed above 10^^ eV indicate a nuclear 
composition (i.e., some mixture of isotopes with mass number A > 4) with 
a relatively small proton component. 

However, the quantitative interpretation of the data in terms of cosmic 
ray composition requires a comparison to shower simulations. Because the 
relevant energies have not been explored in laboratory experiments, the use 
of interaction parameters extrapolated from lower energy experiments is re- 
quired. For instance, a comparison of recent LHC data to predictions of 
different hadronic interaction models showed reasonable overall agreement, 
but each model had its shortcomings |2| . This model uncertainty introduces 
an uncertainty in the interpretation of the X^ax measurements. 

There are also other results that are less suggestive of a nuclear compo- 
sition at the highest energies. The HiRes Collaboration has reported results 
|3| that are compatible with a pure proton beam. In addition, the correlation 
of the arrival directions of the highest energy cosmic rays with nearby extra- 
galactic matter observed by the Auger Observatory j^ is not unexpected if 
the particles are protons, for which magnetic deflections of a few degrees are 
plausible, yet it is more difficult to explain if the particles are heavy nuclei. 

Thus, the current situation is unclear. Additional independent experi- 
mental results related to composition appear mandatory for further progress. 

The muon shower size A*"^ is known to be well correlated to the mass 
of the primary cosmic ray, and a measurement of A^„ would generally be 



independent of the measurement of X^ax- In this way, it is natural to pur- 
sue the use of A^^ in composition studies. However, because of significant 
model uncertainties, obtaining robust composition information from N^^ is 
not straightforward. To briefly demonstrate the problem, consider that the 
difference between the EPOS |5|, |6| and Sibyll 2.1 |7| hadronic interaction 
models with regard to the average muon shower size (A^^) for 10^^ eV proton 
showers is ~ 40%. For comparison, the expected difference between proton 
and iron showers with regard to (A^^) is also ~ 40%. Thus, at present, it is 
quite difficult to make a reliable inference from an observed value of (iV^) to 
the average nuclear mass. 

This does not mean that all useful composition information contained in 
the N^ observable is muddled by model uncertainties, but care must be taken 
as to what inferences to draw. There may be other useful inferences about 
the primary composition besides the average nuclear mass. 

Here we investigate what can be learned from the statistical correlation 
between the two observables Xmax and A^^ assuming a set of events where 
both quantities are measured simultaneously. Our definition of correlation 
(precisely defined in Section 2) contains the lowest-order moment of a two- 
dimensional distribution taking both observables into account, and it is in- 
dependent of systematic uncertainties of the absolute scales of the two ob- 
servables. To our knowledge, such a study has not been described in detail 
in the literature. 

We find that the correlation between X^ax and A"^ can provide a measure 
of the spread of masses in the primary beam. Particularly, one can differ- 
entiate between a well- mixed composition (i.e., a beam with large fractions 
of both light and heavy primaries) and a relatively pure composition (i.e., a 
beam that contains species all of a similar mass). This differentiation remains 
significant and meaningful, even though diluted, in the presence of model un- 
certainties. Moreover, a transition of composition (e.g., from light to heavy) 
will show up as a characteristic change of the correlation with energy. Thus, 
while not providing a measure of the absolute average mass of cosmic rays, 
we find the correlation between X^ax and N^ to be a fairly robust measure 
of significant, complementary characteristics of the cosmic-ray composition. 

The structure of the paper is as follows. In Sec. [2l we define the observ- 
ables and the correlation factor. The analysis in Sec. [3] first regards the case 
of an ideal detector, and then the effects of a realistic detector, of model un- 
certainties and of different composition mixtures. Sec. H] provides a further 
discussion of the method. The paper is concluded in Sec. O 



2. The Xmax and AT^ shower observables 

For our discussion here, we assume that X^ax and A''^ are measured with 
different detection systems, such that the measurement errors are not cor- 
related. For example, X^ax could be measured with an array of telescopes 
that detect the fluorescence and/or Cherenkov light created by the air shower, 
and Nfj_ could be measured by an array of muon detectors deployed at ground 
level. These techniques have been successfully used to measure showers in 
the UHE range 0, i, 0, [ill, [l2 . 



We assume that both measurement techniques are employed on the same 
set of showers. As in Ref. |l[, we define X^ax to be the atmospheric depth 
where the energy deposit is maximum. We define A*"^ to be the total number 
of muons above 1 GeV that reach ground level. Our results are not strongly 
affected by this choice of an energy cut. 

It is typical for the muon shower size to be quoted as the number of 
muons per unit area at a certain distance from the shower axis, with the 



distance being optimized for a particular detector (e.g., [ll|, ll2|). Such a 
muon density observable is highly correlated with A^^. We use Nfj_ because 
it is the more general shower observable (i.e., the specification of a certain 
shower distance is not required). Furthermore, an accurate simulation of 
the lateral development of the shower is not required. This removes any 



potential inaccuracies introduced by thinning algorithms [13|, [Ij] (common 
to many 3-dimensional air shower simulation codes). 

The muon shower size must also be specified at a certain atmospheric 
depth. Here, we assume for simplicity that A*"^ is measured at a constant 
slant depth of 1200 g/cm^, which is well past the range of X^ax observed 
by the Auger Observatory [l|. We checked that the main findings of this 
study do not depend on the specific choice of the observation level of A"^, 
as long as A"^ is observed past shower maximum. We take the zenith angle 
of the showers to be a constant 45°, which corresponds to a slant depth of 
~ 1200 g/cm^ if the muon detector is located ~ 1400 m above sea level. 
In reality, there would be a distribution of zenith angles (or equivalently, of 
slant depths). In this case, the dependence of A"^ on zenith angle should 
be removed by normalizing each A"^ measurement to a central zenith angle 
using a parametrized function derived from the complete A^^ data set. 

We define the correlation factor r for a set of showers to be the linear 



correlation between the observed values of Xmax and A^^ 



COv(Xmax, N^) 



(1) 



where cov and a denote the covariance and standard deviation operators, 
respectively. 

3. Analysis 



For our analysis we use the Conex p^, |16[ air shower simulation package. 
Conex uses a state-of-the-art hybrid calculation scheme 17|]. That is, it uses 
an explicit Monte Carlo simulation for the highest energy interactions in the 
air shower (i.e., the first several interactions) and nuclear-electro-magnetic 
cascade equations for the lower energy interactions. This allows for a fast, 1- 
dimensional air shower simulation. Conex has been shown to reproduce well 



the results [l5| (e.g., the mean and the natural shower-to-shower fluctuation 



of the number of showers particles vs. depth) of the CORSIKA 18|] air shower 
simulation code. 

3.1. Ideal Detector 

Let us first consider the case where the detectors are ideal, i.e., the detec- 
tors have zero measurement uncertainty. In Fig. [T^, we show a scatter plot 
of Xmax and N^ values for air showers initiated by protons and iron nuclei 
with energy of 10^^ eV (1000 events each). The showers were simulated with 
Conex using the QGSJET-01 [l9| high energy interaction model. 

With the consideration of an ideal detector, protons and iron nuclei are 
well separated on the A^-Xmax plane. Proton showers produce the broader 
distribution; iron showers produce the more narrow distribution. For the pro- 
ton distribution, the correlation factor is r^ = 0.0. For the iron distribution, 
rpe = 0.7. For the union of the two sets, Tp+pf. = —0.51. 

The main result we want to elucidate is the following. For composition 
mixtures that contain large fractions of two or more chemical species that 
are well separated in mass number (i.e., significant fractions of both light and 
heavy nuclei - a well- mixed composition), the value of r will be significantly 
negative. This is in contrast to when the composition is comprised of only 
one chemical species, for which case the value of r is near zero or positive. 

Note that the result that r is negative for well-mixed compositions de- 
pends mainly on the relative separation between protons and heavy nuclei 



with regard to the two shower observables. For the same initial conditions, 
-^max is ~ 100 g/cm^ deeper for proton showers than for iron showers, and 
A*"^ is ~ 40% greater for iron showers than for proton showers. Both of these 
expectations are not strongly dependent on the details of the high energy 
hadronic interactions. In this way, the connection between a negative value 
of r and the composition mixture is not strongly model dependent. 

As an aside, we note there is also a "geometrical" correlation between 
Xmax and Nfj_: For measurements past the shower maximum at fixed slant 
depth, Nf^ will increase for deeper Xmax showers. This effect shows up par- 
ticularly well in the case of iron and its relatively large value rpe = 0.7. The 
effect will be diluted, however, for realistic detector conditions. 

3.2. Realistic detector 

We take our nominal Xmax measurement uncertainty to be 20 g/cm^. 
This is similar to the measurement uncertainty of X^ax with the fluorescence 
detector of the Pierre Auger Observatory [1.]. We take our nominal N^ mea- 
surement uncertainty to be 20%. This is similar to the resolution of the 



muon density observable p^ used with the Yakutsk [ll| detector arrays. The 
AMIGA \12^ detector array is expected to have a similar resolution. 

In Fig. [T)d, we show a scatter plot of Xmax and X^ values, which includes 
our nominal measurement uncertainties. As before, the distributions are for 
air showers initiated by protons and iron nuclei with an energy of 10^^ eV. 

In contrast to the ideal detector consideration, protons and iron nuclei 
are not well separated in Fig. [T]d. For the proton distribution, the value of r 
remains Vp = 0.0. For the iron distribution, the value of r decreases, rpe = 
0.1. For the union of the two sets, the value of r increases, r^+^e = —0.32. 

Because X^ax and X^ are measured with independent detectors, the ef- 
fect of measurement uncertainty is to de-correlate the observables. That 
is, the values of r are closer to zero than in the case of the ideal detector. 
However, the main observation from before remains valid: the value of r can 
discriminate between a well-mixed composition and a pure composition. 

The discrimination power of r holds also for data sets smaller than the 
one shown in Fig. 1. For example, for a data set of 200 events, the 1-sigma 
statistical uncertainty on r is ~ 0.07. For this case, r^ and Vp+pe are separated 
by nearly 5-sigma. 



3. 3. Effects of measurement resolution and model uncertainties 

The ability of r to discriminate between composition mixtures is depen- 
dent on the measurement uncertainties. In Fig. |2l we show the value of r for 
different composition scenarios versus the measurement uncertainty of the 
N^ observable, keeping the X^ax uncertainty at 20 g/cm^. The results are 
for 10^^ eV showers. Fig. 2a is calculated with the QGSJET-01 interaction 
model, while Fig. 2b is calculated with the Sibyll 2.1 interaction model. Out 
of the four models considered in this paper - QGSJET-01, QGSJET 11-03 



201, Sibyll 2.1, and EPOS - QGSJET-01 and Sibyll 2.1 give the minimum 
and maximum values of r^+j^e, respectively. 

Fig. 2 demonstrates that the separation between a pure composition and 
a well-mixed composition is significant over a broad range of N^ measurement 
uncertainties. 

Fig. 2 also demonstrates how the value of r depends on the interaction 
model used in the simulation. Comparing QGSJET-01 and Sibyll 2.1, the 
values of r are similar for both helium and iron nuclei. Indeed, we have 
found that for all pure beams with A > 4, r is similar for all four interaction 
models, when realistic measurement uncertainties are taken into account. 

However, the situation for protons is different. The value of r^ is some- 
what model dependent. Related to this, the value of r for mixtures that 
contain a significant proton fraction is also model dependent. The value of 
rp calculated with either the Sibyll 2.1 or EPOS models is significantly more 
negative than the value of rp calculated with either of the QGSJET models. 
Without going into the details of the hadronic interactions themselves, we 
offer the following observation as an insight into the cause of this difference. 
For the Sibyll 2.1 and EPOS models, relative to the QGSJET models, it is 
more common for a large faction of the primary's energy to be transferred to 
electromagnetic particles early on in the development of the shower. This fea- 
ture tends to produce proton showers that are somewhat photon-like (i.e., a 
deep Xjnax and small A^^). In turn, the presence of these photon-like showers 
tends to make the value of r more negative for protons. 

As also discussed below, the model uncertainty reduces the discrimination 
power of r, but useful constraints are still possible. For instance, with our 
nominal measurement uncertainties and for a data set of 200 events, an 
observation of r < —0.25 will favor a well-mixed composition independent of 
interaction models. 



3.4- Sensitivity to different composition mixtures 

In Fig. ini we show the value of r for a proton-iron mixture as a function 
of the number ratio p/{p + Fe). The value of r is calculated with our nominal 
measurement uncertainties. We show results for all four interaction models. 

The minimum value of r occurs for a ratio p/{p + Fe) ^ 0.5. Near this 
ratio the value of r is fairly flat. For example, in the range 0.2 < p/{p+Fe) < 
0.8, r changes by < 0.1. This is true for all models. For well-mixed proton- 
iron mixtures, the value of r can differ by up to ~ 0.18 between the models. 
At higher (lower) proton fractions, r increases sharply and approaches Vp 
{rpe)- A value of r near zero is indicative of a fairly pure composition, while 
a more negative value of r occurs only for a well-mixed beam. 

In Fig. m we show the correlation factor for several different composi- 
tion mixtures calculated with our nominal measurement uncertainties. In 
Fig. H^, the calculations were performed with the QGSJET-01 interaction 
model, while in Fig. |3)d the calculations were performed with the Sibyll 2.1 
interaction model. In most cases, these two models bracket the values of 
r calculated with the QGSJET 11-03 and EPOS models. We plot r as a 
function of RMS(ln(74)) = ■\/var(ln(A)), where A is the mass number of the 
cosmic rays present in the beam and var is the variance operator. The quan- 
tity RMS(ln(74)) is a measure of the purity of the beam (i.e., a pure beam 
will have RMS(ln(A)) = 0, and a mixed beam will have RMS(ln(A)) > 0). 

We show four different pure beams: p. He, C, and Fe; four different bi- 
species beams: p & C, p & Fe, He & Fe, and C & Fe; and one quad-species 
beam: p & He & C & Fe. For each bi-species beam, we plot three different 
mixture ratios: 80%-20%, 50%-50%, and 10%-90%, where the first percentage 
is the fraction of the light component. For the quad-species beam, we test 
the case where all four species are present in equal measure. 

Over a wide variety of mixtures, the correlation between r and RMS(ln(y4)) 
is quite strong: r decreases as RMS(ln(A)) increases. The most negative 
value of r occurs when the mixture is dominated by nearly equal portions of 
proton and iron nuclei. 

From an inspection of Fig. HI we see that an observation of r < —0.25 
would indicate a well-mixed beam, i.e., RMS(ln(y4)) > 1.3. In turn, an ob- 
servation of r > —0.05 would indicate a fairly pure beam, i.e., RMS(ln(74)) < 
0.7. This indication is independent of the hadronic model. 

In contrast, such a model independent indication of the composition seems 
impossible to achieve from the N^ observable alone (i.e., the model uncertain- 
ties dominate everything). An example of this was given in the introduction. 
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Let us give one more example here. We show in Fig. 5 the RMS spread in 
A*"^ for the same composition mixtures shown in Fig. 4. While there might 
be a moderate correlation between RMS(A''^) and RMS(ln(y4)) for a given 
model, the difference in RMS(A^^) between models clearly dominates. 

Given a 2-d distribution of A^^ and Xmax measurements, the lowest order 
moment involving A*"^ that has a clear interpretation regarding the com- 
position seems to be the normalized covariance (i.e., r). There are, of 
course, other ways to look at this problem. One could form a linear com- 
bination like aXinax + hN^. The second moment of this new variable con- 
tains cov(Xniax, X^), and, if formulated carefully, one should expect this new 
variable to have some similarities to r with regard to composition analy- 
sis. Notwithstanding the full range of more complex formalizations, we have 
found the r variable both to be simple and to have a clear, meaningful inter- 
pretation. 

4. Further discussions 

We have investigated the use of the Spearman rank correlation coefficient 
as a substitute for r. We found no significant improvement. 

In Section 3, we only considered showers with an energy of 10^^ eV. 
We checked that the dependence of r on the composition mixture is similar 
throughout the UHE regime. 

When calculating r from data (i.e., not simulations), the showers will be 
distributed within an energy bin. In this case, the X^ax and X^ observables 
must be normalized to a central energy. If the energy estimate for the shower 
has a statistical uncertainty of less than 10% (expected if the longitudinal 



profile of the shower is well measured 2l|), this normalization procedure does 
not add appreciably to the measurement uncertainty of either X^ax oi N^. 

As an observable, the value of r is rather robust against systematic errors 
associated with the detectors. For example, a constant offset or multiplicative 
factor in the measurement of X^ax or X^ does not affect the calculation of 
r. 

If there is a light to heavy composition transition, as indicated by recent 
elongation rate measurements, then there will likely be an energy range where 
the composition is well-mixed. The work presented here indicates that it is 
possible to obtain a robust, independent indication of a well- mixed condition. 
The discrimination power of r is best between a pure composition and a 
composition with significant fractions of both protons and iron nuclei. It is 



interesting that such a well-mixed scenario, near 10 eV, has been suggested 



in the literature 22 



It is possible to measure r as a function of energy (i.e., in consecutive 
energy bins). If there exists good statistics over perhaps an order of magni- 
tude in energy, it may be possible to observe the composition entering and/or 
exiting a well-mixed state, i.e., to observe r changing with energy. As an ex- 
ample, consider the case where the composition is transitioning from a nearly 
pure beam of protons (or iron) to a well-mixed beam of protons and iron. In 
this case, the expected change of r with energy can be obtained from Fig. [31 
The observation of such a change in r with energy would be a particularly 
strong, model independent indication of a composition transition. 

For a given interaction model to provide a successful description of the 
shower data, there must be consistency between the range of RMS(ln(y4)) im- 
plied by the observed (Xmax) and RMS(Xinax) and the range of RMS(ln(y4)) 
implied by the observed r. For example, suppose that (Xmax) is observed to 
be somewhat less than (-'^max)„roton ^^^ ^^^^ RMS(Xinax) is observed to be 
broad such that the only composition mixtures that fit well are well-mixed. 
Then if the observed value of r is inconsistent with a well-mixed composi- 
tion, the interaction model cannot be judged to be self-consistent. In this 
way, a measurement of r can be used to evaluate the self-consistency of the 
interaction models. 

An observed value of r inconsistent with a pure composition, in the frame- 
work of the models discussed here, does not exclude the possibility that the 
composition is actually quite pure and that there are physical processes that 
are currently not accounted for in air shower simulations. However, the 
observed value of r would put constraints on any hypothesized physical pro- 
cesses that are proposed to explain, e.g., the observed elongation rate. 

5. Conclusions 

We have studied in detail the sensitivity of the correlation factor r be- 
tween Xmax and A^^ to qualities of the UHE cosmic ray composition. The 
correlation factor provides composition information that complements the 
standard composition analysis based on the average and RMS of the observ- 
ables. It incorporates the muon shower size observable A^^ and thus provides 
an indication of the composition that is experimentally independent of the 
elongation rate analysis. While not providing a direct indication of the aver- 
age nuclear mass, the value of r is sensitive to whether the UHE cosmic ray 
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beam is pure or well-mixed. Testing whether the beam is pure or well-mixed 
is well motivated by recent measurements of the depth of shower maximum. 
The discrimination power of r is affected by model uncertainties; however, it 
is possible to make useful constraints independent of the interaction models. 
More sophisticated analysis methods might be developed in the future that 
make a combined use of the information contained in the moments of the two- 
dimensional Xjnax - Nf^ - distribution, particularly if the detector effects of a 
specific experimental setup are well understood and if the theoretical model 
uncertainties can be reduced. Here we showed what kind of new information 
is contained in the correlation parameter r. 
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Figure 1: Xmax - Np, distribution for proton showers and iron showers with an energy of 
10^^ eV and a zenith angle of 45°. There are 1000 iron showers and 1000 proton showers 
plotted. The showers were simulated with Conex using the QGSJET-01 high energy 
interaction model. Figure (a) is for an ideal detector (i.e., zero measurement uncertainty). 
Figure (b) is for a realistic detector; see the text for a description of the measurement 
uncertainties. 
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Figure 2: Correlation factor r versus the N^ measurement uncertainty for three different 
pure compositions and one well-mixed composition, a proton-iron mixture. The proton- 
iron mixture is 50% protons and 50% iron nuclei. The results in Figure (a) are derived 
with the QGSJET-01 interaction model. The results in Figure (b) are derived with the 
Sibyll 2.1 interaction model. 
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Figure 3: Correlation factor r for a proton-iron composition mixture. The quantity p/(j) + 
Fe) is the ratio of the number of protons in the beam to the total number of particles. 
We show the results of four different interaction models. We assume a realistic detector 
with the nominal measurement uncertainties described in the text. 
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Figure 4: Correlation factor r for several different composition mixtures versus 
RMS(ln(A)). For each bi-species beam, we plot three different mixture ratios: 80%-20%, 
50%-50%, and 10%-90%, where the first percentage is the fraction of the light component. 
For the quad-species beam, we test the case where all four species are present in equal 
measure. We assumed a realistic detector with the nominal measurement uncertainties 
described in the text. 
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Figure 5: RMS(A''^) for several different composition mixtures versus RMS(ln(A)). We 
assumed a realistic detector with the nominal measurement uncertainties described in the 
text. 
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